### Preliminary ----------------------------------------------------------------
# Load libraries
library(cowplot)
library(fixest)
library(janitor)
library(tidyverse)
library(magrittr)
library(bayestestR)
library(lubridate)

# Load data
slips_uncleaned <- read_csv("Slip Data (Original Download).csv")

load("Bill Dataframe.Rda")
load("Ideal Points Dataframe.Rda")
load("Witness Slips Dataframe (Cleaned).Rda")
load("Witness Dataframe.Rda")

ideal_sm <- read_tsv("sm ideal points.tab")

### Figure 1: Distribution of Slips per Bill and Support -----------------------
## Number of Slips Per Bill
# Format data for graphing
slips_per_bill <- slips %>%
  group_by(bill_name, session) %>%
  summarize(n_witnesses = n())

# Plot distribution of the number of witnesses per bill
volume_distribution <- ggplot(slips_per_bill, aes(x = n_witnesses)) +
  geom_histogram(bins = 20, fill = "lightgrey", color = "black") +
  geom_vline(aes(xintercept = median(slips_per_bill$n_witnesses)),
             lty = "dashed") + 
  scale_x_log10(breaks = scales::log_breaks(n = 5),
                labels = scales::comma, 
                expand = c(0.01, 0.01)) +
  scale_y_continuous(expand = c(0,0), limits = c(0, 4500)) + 
  labs(x = "Witness Slips per Bill (Log Scale)", y = "Count") + 
  theme_bw()  +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank()
  )

## Distribution of Support by Bill
# Format data for graphing
support_df <- slips_uncleaned %>%
  summarize(Favor = sum(position == "prop") / n(),
            Oppose = sum(position == "opp") / n(),
            `No Position` = sum(position == "nopos") / n()) %>%
  pivot_longer(cols = c(Favor, Oppose, `No Position`),
               names_to = "Position", 
               values_to = "Percent")

# Plot distribution of support
support_distribution <- ggplot(support_df, aes(x = Position, fill = Position, y = Percent)) + 
  geom_col(position = position_dodge(), color = "black") +
  scale_y_continuous(labels = scales::percent, expand = c(0,0), limits = c(0, 0.65)) + 
  scale_fill_manual(values = c("white", "gray", "black")) + 
  labs(x = "Position", y = "Percent of All Slips", fill = NULL) +
  theme_bw() +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.grid.major.x = element_blank(),
    legend.position = "none"
  )

# Combine plots
plot_grid(volume_distribution, support_distribution)

### Figure 2: Partisan Direction of Witness Slip Support ------------------------
## Partisanship of support (bill-level)
# Calculate net support by party
bill_partisanship <- slips %>%
  group_by(session, bill_name, sponsor_party) %>%
  summarize(support = sum(position == 1, na.rm = TRUE),
            opposition = sum(position == -1), na.rm = TRUE,
            net_support = sum(position, na.rm = TRUE)) %>% 
  group_by(sponsor_party) %>%
  summarize(av_support = mean(support, na.rm = TRUE),
            av_oppose = mean(opposition, na.rm = TRUE),
            av_net_support = mean(net_support, na.rm = TRUE),
            se_net_support = sd(net_support, na.rm = TRUE) / sqrt(n())) %>%
  filter(sponsor_party %in% c("D", "R"))

# Put data into form for graphing
graphing_data <- data.frame(
  coefficient = c("Democratic Bills", "Republican Bills"),
  estimate = bill_partisanship$av_net_support,
  se = bill_partisanship$se_net_support) %>%
  mutate(lower = estimate - 1.96 * se,
         upper = estimate + 1.96 * se)

# Create bill-level plot
bill_level <- ggplot(graphing_data, aes(x = estimate, y = coefficient)) +
  geom_pointrange(aes(xmin = lower, xmax = upper)) +
  geom_vline(aes(xintercept = 0), lty = "dashed") +
  #scale_x_continuous(limits = c(-23, 23)) +
  labs(x = "Net Support \n (Supporting Slips - Opposing Slips)", 
       y = NULL, title = "Bill-Level Partisanship") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5))

## Partisanship of support (witness-level)
# Put data into form for graphing
graphing_data <- data.frame(
  coefficient = "    All Witnesses",
  estimate = mean(witnesses$net_support, na.rm = TRUE),
  se = sd(witnesses$net_support, na.rm = TRUE) / 
    sqrt(nrow(witnesses))) %>%
  mutate(lower = estimate - 1.96 * se,
         upper = estimate + 1.96 * se)

# Create plot
witness_level <- ggplot(graphing_data, aes(x = estimate, y = coefficient)) +
  geom_pointrange(aes(xmin = lower, xmax = upper)) +
  geom_vline(aes(xintercept = 0), lty = "dashed") +
  scale_x_continuous(limits = c(-0.5, 0.5)) +
  labs(x = "Net Partisanship \n (Republican Positions - Democratic Positions)", 
       y = NULL, title = "Witness-Level Partisanship") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5))

# Combine plots
plot_grid(bill_level, witness_level, nrow = 2)

### Figure 3: Ideology of Witnesses and Slips, Relative to State Legislators----
# Prepare Shor-McCarty data
ideal_sm <- ideal_sm %>%
  select(name, party, st, np_score, 
         senate2012:senate2020,
         house2012:house2020) %>%
  mutate_at(vars(senate2012:house2020), function(x) ifelse(is.na(x), 0, x)) %>%
  filter(if_any(senate2012:house2020) > 0)

# Calculate party medians for IL legislators and other quantities
median_dem_ill <- ideol_estimates %>%
  filter(party == "D") %>%
  pull(D1) %>%
  median(na.rm = TRUE)

median_rep_ill <- ideol_estimates %>%
  filter(party == "R") %>%
  pull(D1) %>%
  median(na.rm = TRUE)

min_dem_ill <- ideol_estimates %>%
  filter(party == "D") %>%
  pull(D1) %>%
  min(na.rm = TRUE)

max_rep_ill <- ideol_estimates %>%
  filter(party == "R") %>%
  pull(D1) %>%
  max(na.rm = TRUE)

median_leg_ill <- ideol_estimates %>%
  pull(D1) %>%
  median(na.rm = TRUE)

# Calculate party medians for all state legislators
median_dem_all <- ideal_sm %>%
  filter(party == "D") %>%
  pull(np_score) %>%
  median(na.rm = TRUE)

median_rep_all <- ideal_sm %>%
  filter(party == "R") %>%
  pull(np_score) %>%
  median(na.rm = TRUE)

median_all <- ideal_sm %>%
  pull(np_score) %>%
  median(na.rm = TRUE)

# Calculate min and max for IL legislators
min_dem_ill <- ideol_estimates %>%
  filter(party == "D") %>%
  pull(D1) %>%
  min(na.rm = TRUE)

max_rep_ill <- ideol_estimates %>%
  filter(party == "R") %>%
  pull(D1) %>%
  max(na.rm = TRUE)

# Calculate party midpoint
party_midpoint_ill <- (median_dem_ill + median_rep_ill) / 2
party_midpoint_all <- (median_dem_all + median_rep_all) / 2

# Calculate % witnesses located closer to Republicans versus Democrats
(ideol_estimates$D1[ideol_estimates$n_slips >= 20] > party_midpoint_ill) %>%
  mean(na.rm = TRUE)

(ideol_estimates$D1[ideol_estimates$n_slips >= 20] > party_midpoint_all) %>%
  mean(na.rm = TRUE)

(ideol_estimates$D1[ideol_estimates$n_slips >= 20] < min_dem_ill) %>%
  mean(na.rm = TRUE)

(ideol_estimates$D1[ideol_estimates$n_slips >= 20] > max_rep_ill) %>%
  mean(na.rm = TRUE)

# Restrict data to witnesses and rename variables
witness_df <- ideol_estimates %>%
  filter(n_slips >= 20,
         type == "Witnesses") %>%
  mutate(category = case_when(
    group_affil == "Group-Affiliated" ~ "Group-Affiliated",
    group_affil == "Non-Affiliated" ~ "Non-Affiliated"),
    category = factor(category, levels = c("Group-Affiliated",
                                           "Non-Affiliated")))

# Add rows to weight by the number of slips
witness_df <- witness_df %>% 
  filter(n_slips >= 20,
         type == "Witnesses") %>%
  uncount(n_slips) %>%
  mutate(type = "Witnesses (Weighted)") %>%
  bind_rows(witness_df, .)

# Used regression models to estimate means with clustered standard errors
model_1 <- feols(D1 ~ 1,
                 cluster = "ids",
                 data = filter(witness_df, type == "Witnesses"))

model_2 <- feols(D1 ~ 1,
                 cluster = "ids",
                 data = filter(witness_df, type == "Witnesses (Weighted)"))

# Put estimates into dataframe for graphing
graphing_data <- data.frame(
  type = c("Witnesses", "Slips"),
  Mean = c(coef(model_1), coef(model_2)),
  SE = c(se(model_1), se(model_2))) %>%
  mutate(Lower = Mean - 1.96 * SE,
         Upper = Mean + 1.96 * SE)
  
# Create plot
graphing_data %>%
  ggplot(aes(x = Mean, y = type)) +
  geom_pointrange(aes(xmin = Lower, xmax = Upper), size = 0.2) +
  geom_text(aes(label = round(Mean, 2)), vjust = -0.5, size = 3) +
  scale_x_continuous(limits = c(median_dem_ill - 0.1, median_rep_ill + 0.1),
                     breaks = c(median_dem_ill, party_midpoint_ill, 0, median_rep_ill),
                     labels = c("Median Democrat \n (-0.76)", "Party Midpoint \n (-0.13)", "\n (0)", "Median Republican \n (0.50)")) +
                     #labels = c("Median \n Democrat", "Party \n Midpoint", "Median \n Republican")) + 
  scale_shape_manual(values = c(21, 16),
                     guide = guide_legend(reverse = TRUE)) + 
  labs(x = "Shor-McCarty (NP) Ideal Point", y = NULL, shape = NULL) + 
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    legend.position = "bottom"
  )

### Figure 4: Ideological Distributions of Legislators and Witnesses -----------
# Create density plot
graphing_data <- ideol_estimates %>%
  filter(n_slips >= 20 | type == "Legislators") %>%
  mutate(party = ifelse(party == "D", "Democrats", "Republicans"),
         group = ifelse(type == "Legislators", party, type))

ggplot(graphing_data, aes(x = D1, fill = group, lty = group)) +
  geom_density(alpha = 0.2, adjust = 1.5) +
  facet_wrap(~type, scales = "free_y", nrow = 3) + 
  scale_fill_manual(values = c("skyblue", "lightcoral",  "yellow4")) +
  scale_linetype_manual(values = c("solid", "dashed", "solid")) +
  scale_x_continuous(breaks = seq(-4, 4, 2),
                     limits = c(-4.2, 4.2)) + 
  scale_y_continuous(expand = c(0,0)) + 
  labs(x = "Ideal Point Estimates", y = "Density",
       fill = NULL, linetype = NULL) + 
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    legend.position = "bottom"
  )

### Figure 5: Differences in Group- versus Non-Affiliated Witnesses-------------
## Partisanship of support (bill-level)
# Calculate net support by party
bill_partisanship <- slips %>%
  group_by(session, bill_name, sponsor_party, group_affiliated) %>%
  summarize(net_support = sum(position, na.rm = TRUE)) %>% 
  group_by(sponsor_party, group_affiliated) %>%
  summarize(av_net_support = mean(net_support, na.rm = TRUE),
            se_net_support = sd(net_support, na.rm = TRUE) / sqrt(n())) %>%
  filter(sponsor_party %in% c("D", "R"))

# Put data into form for graphing
graphing_data <- data.frame(
  coefficient = c("Democratic Bills", "Democratic Bills", 
                  "Republican Bills", "Republican Bills"),
  group = bill_partisanship$group_affiliated,
  estimate = bill_partisanship$av_net_support,
  se = bill_partisanship$se_net_support) %>%
  mutate(lower = estimate - 1.96 * se,
         upper = estimate + 1.96 * se)

# Create bill-level plot
bill_level <- ggplot(graphing_data, aes(x = estimate, y = coefficient)) +
  geom_pointrange(aes(xmin = lower, xmax = upper, color = group, shape = group),
                  position = position_dodge(width = 0.5)) +
  geom_vline(aes(xintercept = 0), lty = "dashed") +
  #scale_x_continuous(limits = c(-23, 23)) +
  scale_color_manual(values = c("black", "gray")) +
  labs(x = "Net Support \n (Supporting Slips - Opposing Slips)", 
       y = NULL, title = "a) Bill-Level Partisanship",
       color = NULL,
       shape = NULL) + 
  theme_bw() +
  theme(
    legend.position = "none",
    plot.title = element_text(hjust = 0.5)
  )

## Partisanship of support (witness-level)
# Put data into form for graphing
graphing_data <- witnesses %>%
  group_by(group_affiliated) %>%
  summarize(estimate = mean(net_support, na.rm = TRUE),
            se = sd(net_support, na.rm = TRUE) / sqrt(n()))  %>%
  mutate(lower = estimate - 1.96 * se,
         upper = estimate + 1.96 * se)

# Create plot
witness_level <- ggplot(graphing_data, aes(x = estimate, y = group_affiliated)) +
  geom_pointrange(aes(xmin = lower, xmax = upper, shape = group_affiliated, color = group_affiliated)) +
  geom_vline(aes(xintercept = 0), lty = "dashed") +
  scale_color_manual(values = c("black", "gray")) +
  scale_x_continuous(limits = c(-1, 1)) +
  labs(x = "Net Partisanship \n (Republican Positions - Democratic Positions)", 
       y = NULL, title = "b) Witness-Level Partisanship",
       color = NULL,
       shape = NULL) + 
  theme_bw() +
  theme(
    legend.position = "none",
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank(),
    plot.title = element_text(hjust = 0.5)
  )

## Witness and slip ideology plot
# Estimate models
model_3 <- feols(D1 ~ -1 + group_affil,
                 cluster = "ids",
                 data = filter(witness_df, type == "Witnesses"))

model_4 <- feols(D1 ~ -1 + group_affil,
                 cluster = "ids",
                 data = filter(witness_df, type == "Witnesses (Weighted)"))

# Put data into form for graphing
graphing_data <- data.frame(
  type = c("Witnesses", "Witnesses", "Slips", "Slips"),
  affiliation = c("Group-Affiliated", "Non-Affiliated", 
                  "Group-Affiliated", "Non-Affiliated"),
  Mean = c(coef(model_3), coef(model_4)),
  SE = c(se(model_3), se(model_4))) %>%
  mutate(Lower = Mean - 1.96 * SE,
         Upper = Mean + 1.96 * SE)

# Create main plot
ideol_plot <- graphing_data %>%
  ggplot(aes(x = Mean, y = type, color = affiliation, shape = affiliation)) +
  geom_pointrange(aes(xmin = Lower, xmax = Upper)) +
  scale_x_continuous(limits = c(median_dem_ill - 0.1, median_rep_ill + 0.1),
                     breaks = c(median_dem_ill, party_midpoint_ill, 0, median_rep_ill),
                     labels = c("Median Dem. \n (-0.76)", "Party Midpoint \n (-0.13)", "\n (0)", "Median Rep. \n (0.50)")) +
  scale_color_manual(values = c("black", "gray")) +
  labs(x = "Shor-McCarty (NP) Ideal Point", y = NULL, shape = NULL, color = NULL, 
       title = "c) Witness and Slip Ideology") + 
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  )

## Distribution of witness ideology
# Prepare data for graphing
graphing_data <- ideol_estimates %>%
  filter(n_slips >= 20) %>%
  filter(type == "Witnesses") %>%
  mutate(category =  group_affil,
         facet = "Witnesses")

# Create density plots
ideol_dist <- ggplot(graphing_data, aes(x = D1, fill = category, linetype = category)) +
  geom_density(alpha = 0.2, adjust = 1.5) +
  scale_fill_manual(values = c("black", "lightgray")) +
  scale_linetype_manual(values = c(1, 2)) + 
  scale_x_continuous(breaks = seq(-4, 4, 2),
                     limits = c(-4.2, 4.2)) + 
  scale_y_continuous(expand = c(0,0)) + 
  labs(x = "Ideal Point Estimates", y = "Density",
       fill = NULL, linetype = NULL, title = "d) Distribution of Witness Ideology") + 
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  )

# Combine plots into single figure
plot_grid(bill_level,  witness_level,ideol_plot, ideol_dist,
          nrow = 2)

### Figure 6: Witness Slip Involvement and Ideology, Interest Groups -----------
# Create dataframe of ideology by interest group
ig_df <- ideol_estimates %>%
  filter(type == "Witnesses", 
         group_idname != "individual/missing") %>%
  group_by(group_idname) %>%
  summarize(mean_ideol = mean(D1),
            n_witnesses = n(),
            n_slips = sum(n_slips)) %>%
  filter(n_witnesses >= 5) %>%
  ungroup()

# Estimate the percent closer to Democrats than Republicans
(ig_df$mean_ideol > party_midpoint_ill) %>%
  mean(na.rm = TRUE)

# Create dataframe for most active interest groups 
top_igs <- ideol_estimates %>%
  filter(type == "Witnesses", group_idname != "individual/missing") %>%
  group_by(group_idname) %>%
  mutate(mean_ideol = mean(D1),
         n_witnesses = n(),
         n_slips = sum(n_slips),
         group_idname = str_to_title(group_idname),
         group_idname = str_replace(group_idname, "Illinois", "IL"),
         group_idname = str_replace(group_idname, "Aclu", "ACLU"),
         group_idname = str_replace(group_idname, "Unitarian Universalist Advocacy Network Of IL", "Unitarian Universalist Advocacy Network"),
         group_idname = str_replace(group_idname, "Igor", "Illinois Gun Owners' Rights"),
         group_idname = str_replace(group_idname, "Afscme", "State, County and Municipal Employees (AFSCME)"),
         group_idname = str_replace(group_idname, "Education", "Teachers' Unions"),
         group_idname = str_replace(group_idname, "One Northside", "Organizing Neighborhoods for Equality"),
         group_idname = str_replace(group_idname, "Ic4ic", "IL Coalition for Informed Consent")) %>%
  filter(n_witnesses >= 5, n_slips >= 2800) %>%
  ungroup()

# Create plot of average ideology by interest group
ggplot(ig_df, aes(x = mean_ideol, y = n_slips)) +
  geom_point(aes(size = n_witnesses), shape = 1, alpha = 0.5) +
  geom_text(data = top_igs,
            aes(label = group_idname, y = n_slips + 400),
            size = 2.5) +
  labs(x = "Ideal Point Average", 
       y = "Number of Witness Slips \n Filed by Group Members",
       size = "Number of Unique Witnesses") + 
  scale_y_continuous(expand = c(0,0),
                     limits = c(0, 14500),
                     breaks = seq(2500, 12500, 5000)) +
  scale_x_continuous(breaks = seq(-2, 2, 1),
                     limits = c(-2.3, 2.3)) +
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    legend.position = "bottom"
  )

### Figure 7: Distribution of Witness Ideology Within Groups -------------------
# Create dataframe for most active interest groups 
top_igs <- ideol_estimates %>%
  filter(type == "Witnesses", group_idname != "individual/missing") %>%
  group_by(group_idname) %>%
  mutate(mean_ideol = mean(D1),
         n_witnesses = n(),
         n_slips = sum(n_slips),
         group_idname = str_to_title(group_idname),
         group_idname = str_replace(group_idname, "Illinois", "IL"),
         group_idname = str_replace(group_idname, "Aclu", "ACLU"),
         group_idname = str_replace(group_idname, "Igor", "IL Gun Owners' Rights"),
         group_idname = str_replace(group_idname, "Afscme", "State, County and Municipal Employees (AFSCME)"),
         group_idname = str_replace(group_idname, "Education", "Teachers' Unions"),
         group_idname = str_replace(group_idname, "One Northside", "Organizing Neighborhoods for Equality"),
         group_idname = str_replace(group_idname, "Ic4ic", "IL Coalition for Informed Consent"),
         group_idname = str_replace(group_idname, "Dna", "DNA")) %>%
  filter(n_witnesses >= 15) %>%
  ungroup()

# Create plot of witness ideology for most active interest groups
top_igs %>%
  mutate(group_idname = as.factor(group_idname),
         group_idname = reorder(group_idname, D1)) %>%
  ggplot(aes(x = D1)) +
  geom_dotplot(fill = "black", binwidth = 0.05) + 
  facet_grid(group_idname ~ ., scales = "free_y") + 
  labs(x = "Ideal Point", y = "Count") + 
  theme_bw() + 
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    strip.text.y = element_text(angle=0)
  )

### Figure 8: Witness Ideology by Committee Jurisdiction -------------------
# classify committees by jurisdiction
slips <- slips %>%
  mutate(committee_recode = case_when(
      committee == "Executive" ~ "Executive",
      committee == "Judiciary" ~ "Criminal Justice and Judiciary",
      committee == "Human Services" ~ "Health and Human Services",
      committee == "Revenue" ~ "Revenue and Finance",
      committee == "Criminal Law" ~ "Criminal Justice and Judiciary",
      committee == "Education" ~ "Education",
      committee == "Licensed Activities and Pensions" ~ "Licensing and Pensions",
      committee == "Transportation" ~ "Transportation",
      committee == "State Government" ~ "State and Local Government",
      committee == "Local Government" ~ "State and Local Government",
      committee == "Insurance" ~ "Insurance",
      committee == "Higher Education" ~ "Education",
      committee == "Public Health" ~ "Health and Human Services",
      committee == "Revenue & Finance" ~ "Revenue and Finance",
      committee == "State Government & Veterans Affairs" ~ "State and Local Government",
      committee == "Judiciary - Criminal" ~ "Criminal Justice and Judiciary",
      committee == "Environment and Conservation" ~ "Environment and Conservation",
      committee == "Judiciary - Civil" ~ "Criminal Justice and Judiciary",
      committee == "Financial Institutions" ~ "Revenue and Finance",
      committee == "Agriculture" ~ "Agriculture",
      committee == "Transportation: Vehicles & Safety" ~ "Transportation",
      committee == "State Government Administration" ~ "State and Local Government",
      committee == "Health Care Licenses" ~ "Health and Human Services",
      committee == "Licensed Activities" ~ "Licensing and Pensions",
      committee == "Personnel & Pensions" ~ "Licensing and Pensions",
      committee == "Energy and Public Utilities" ~ "Energy and Public Utilities",
      committee == "Commerce and Economic Development" ~ "Labor and Commerce",
      committee == "Health" ~ "Health and Human Services",
      committee == "Labor" ~ "Labor and Commerce",
      committee == "Environment" ~ "Environment and Conservation",
      committee == "Counties & Townships" ~ "State and Local Government",
      committee == "Agriculture & Conservation" ~ "Agriculture",
      committee == "Labor & Commerce Committee" ~ "Labor and Commerce",
      committee == "Elem Sec Ed: School Curric Policies" ~ "Education",
      committee == "Appropriations-Human Services" ~ "Revenue and Finance",
      committee == "Labor & Commerce" ~ "Labor and Commerce",
      committee == "Pensions" ~ "Licensing and Pensions",
      committee == "Agriculture and Conservation" ~ "Agriculture",
      committee == "Labor & Commerce Committee" ~ "Labor and Commerce",
      committee == "Elem Sec Ed: School Curric Policies" ~ "Education",
      committee == "Labor & Commerce" ~ "Labor and Commerce",
      committee == "Pensions" ~ "Licensing and Pensions",
      committee == "Agriculture and Conservation" ~ "Agriculture",
      committee == "Public Utilities" ~ "Energy and Public Utilities",
      committee == "Veterans Affairs" ~ "State and Local Government",
      committee == "Government Accountability/Pensions" ~ "State and Local Government",
      committee == "Business Occupational Licenses" ~ "Licensing and Pensions",
      committee == "Elem Sec Ed: Adm., Lic. & Charter" ~ "Education",
      committee == "Energy & Environment" ~ "Energy and Public Utilities",
      committee == "Cities & Villages" ~ "State and Local Government",
      committee == "Energy" ~ "Energy and Public Utilities",
      committee == "Elem Sec Ed: Licensing, Admin." ~ "Education",
      committee == "Behavioral and Mental Health" ~ "Health and Human Services",
      committee == "Elementary & Secondary Education" ~ "Education",
      committee == "Gaming" ~ "Licensing and Pensions",
      committee == "Government Reform" ~ "State and Local Government",
      str_detect(committee, "Firearm") ~ "Firearms",
      TRUE ~ NA))

by_committee <- slips %>%
  select(bill_name, session, witness_id,committee_recode) %>%
  mutate(bill = paste(ifelse(grepl("House Bill",bill_name),"HB",
                             ifelse(grepl("Senate Bill",bill_name),"SB",NA)),
                      as.numeric(grepl("[[:digit:]]+",bill_name)),sep=""))

by_committee <- ideol_estimates %>%
  filter(n_slips >= 20) %>%
  select(ids, D1) %>%
  rename(witness_id = ids) %>%
  inner_join(by_committee, ., by = "witness_id")

# Used regression models to estimate means with clustered standard errors
model_1 <- feols(D1 ~ as.factor(committee_recode) - 1,
                 cluster = "witness_id",
                 data = by_committee)

# Construct list of committees by number of slips
committee_list <- slips %>%
  group_by(committee_recode) %>%
  filter(!(is.na(committee_recode))) %>%
  summarize(numslips = n()) 

# Put estimates into dataframe for graphing
graphing_data <- data.frame(
  Mean = coef(model_1),
  SE = se(model_1)) %>%
  mutate(Lower = Mean - 1.96 * SE,
         Upper = Mean + 1.96 * SE,
         committee = committee_list$committee_recode,
         numslips = committee_list$numslips) %>%
  arrange(desc(committee_list$numslips)) %>%
  mutate(committee = factor(committee),
         committee = reorder(committee, numslips))

# Calculate committee median ideal points
committee_medians <- ideol_estimates %>%
  select(D1,starts_with("com_")) %>%
  pivot_longer(
    cols = starts_with("com_"),
    names_to = "committee",
    values_to = "value"
  ) %>%
  filter(value == 1) %>%
  group_by(committee) %>%
  summarise(median_ideal_point = median(D1, na.rm = TRUE)) %>%
  mutate(committee = committee_list$committee_recode) %>%
  arrange(desc(committee_list$numslips))

# 
graphing_data$committee_median <- committee_medians$median_ideal_point[
  match(graphing_data$committee,committee_medians$committee)]

# Create plot of witness ideology by committee
ggplot(graphing_data, aes(y = Mean, x = committee)) +
  geom_point(size = 1) + 
  geom_point(aes(y=committee_median,x=committee),size=1,
             pch=4) + 
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2) +
  coord_flip() +
  labs(y = "Mean Witness \n Ideal Point",
       x = NULL) + 
  scale_y_continuous(limits = c(median_dem_ill - 0.3, 1.1),
                     breaks = c(median_dem_ill, party_midpoint_ill,0, median_rep_ill),
                     labels = c("Median \n Democrat", "Party \n Midpoint","0", "Median \n Republican")) + 
  theme_bw() +
  theme(
    panel.border = element_blank(),
    axis.line.x.bottom = element_line(color = "black"),
    axis.line.y.left = element_line(color = "black"),
    axis.ticks.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
  )


